package algorithm.fft;

import java.util.Random;

import edu.emory.mathcs.jtransforms.fft.FloatFFT_1D;

public class DepthFirstFFT {
	private static class Control{
		public int step_child;
		public int size_child;
		public float[] real; 
		public float[] image;
		public int k; 
	}
	//trigonometric function
	//private Control cntrl=new Control();
	int dim; 
	private IFFTNodeFactory f=null; 
	private int order;
	private Trigonometry t; 
	private WorkMemory wm; 
	private int window; 
	public DepthFirstFFT(int dim, int order,int w){
		this.dim=dim; 
		this.window=w; 
		this.order=order; 
		t=new Trigonometry(dim,order);
		wm=new WorkMemory(dim);
		f=new FastFFTNodeFactory(log2(dim)*2+1,dim,t,wm);
		BinaryFFTCache.buildCache(); 
		
	}
	/*
	public static int[] doBitReversal(final int order){
		int N=1<<order; 
		int half_order=(int)Math.ceil(order/2.0); 
		for (int i=0; i<N; i++){
			for (int k=0; k<half_order; k++){
				
			}
		}
	}
	*/
	public void fft2(boolean[] x,float[] real,float[] image){
		Control cntrl=new Control(); 
		cntrl.real=new float[dim];
		cntrl.image=new float[dim]; 
		doFFT(x,real,image,0,1,x.length,x.length, order); 
	}
	private void doFFT(boolean[] x,
			float[] real, float[] image,
			int from, int step, int to,int size,int order){
		if (size==FFTNode.MIN_SIZE){
			return; 
		}		
		int step_child=step<<1; 
		int size_child=size>>1;
		doFFT(x,real,image,from,step_child,to-step,size_child,order-1); 
		doFFT(x,real,image,from+step,step_child,to,size_child,order-1);
		float right_re;  
		float right_im;
		float left_re;  
		float left_im; 
        float wk_re;  
        float wk_im; 
        switch(size){
        case 1: 
        	return;
        case 2:        	 
        	left_re=x[from]?1:0; 
        	right_re=x[from+step]?1:0;
        	real[from]=left_re+right_re;  
        	image[from]=0; 
        	real[from+step]=left_re-right_re;
        	image[from+step]=0;        	
        	return ;
        default:
            int offset_leftChild=from;
        	int dual; 
        	int offset_rightChild=from+step; 
        	int half_size=size>>1;
			int window=Math.min(half_size, 16);
    		for (int k = 0; k<window; k++) {
    			Complex e_kN=t.get(k, size, order);

    			right_re=real[offset_rightChild]; 
    			right_im=image[offset_rightChild]; 
    			    			
    			wk_re=(float)(right_re*e_kN.re-right_im*e_kN.im); 
    			wk_im=(float)(right_re*e_kN.im+right_im*e_kN.re);
    			    			
    			left_re=real[offset_leftChild]; 
    			left_im=image[offset_leftChild]; 
    			wm.real[k]=left_re+wk_re;
    			wm.image[k]=left_im+wk_im;
    			
    			dual=k+half_size;
    			wm.real[dual]=left_re-wk_re;
    			wm.image[dual]=left_im-wk_im;
    			offset_leftChild+=step_child;
    			offset_rightChild+=step_child;
    			 
    		}
    		
    		int m=half_size*step;
    		int offset=from;
    		int pos; 
			for (int k = 0; k<window; k++) {
				
				real[offset]=wm.real[k];
				image[offset]=wm.image[k];
				dual=offset+m; 
				pos=k+half_size; 
				real[dual]=wm.real[pos]; 
				image[dual]=wm.image[pos]; 
				offset+=step;
				
			}		
        }
	}
	public static int log2(int n){
		return (int)(Math.log(n)/Math.log(2))+1; 
	}
	public final void fft(boolean[] x, float[] real,float[] image){ //suppose x.length is pow of 2; 
		FFTNode stackTop=f.createRootNode(0, x.length,order);
		stackTop.prev=null; 
		FFTNode[] childNodes=new FFTNode[2];
		int numOfChildren=0; 
		int i; 
		while (stackTop!=null){			   
			if (stackTop.isReady){
				stackTop.compute(x,real,image,window); 
				f.freeNode(stackTop);
				stackTop=stackTop.prev; 
			}else {
				if (stackTop.size>FFTNode.MIN_SIZE){
					numOfChildren=stackTop.createChildNodes(childNodes);
					for (i=0; i<numOfChildren; i++){
						childNodes[i].prev=stackTop; 
						stackTop=childNodes[i];   
					}
				}
			}
		}
	}
	public FFTNode fft3(boolean[] x, float[] real,float[] image){ //suppose x.length is pow of 2; 
		int size=order*2+1;  
		FFTNode[] stack=new FFTNode[size]; // this size should be big than log2(Nlog2N)+1
		int top=0; 
		//float[] real=new float[x.length]; 
		//float[] image=new float[x.length]; 
		stack[top++]=f.createRootNode(0, x.length,order);
		FFTNode[] childNodes=new FFTNode[2];
		int numOfChildren=0; 
		while (top>0){
			FFTNode current=stack[top-1];  
			if (current.isReady){
				current.compute(x,real,image,window); 
				f.freeNode(current); 
				top--; 
			}else {
				if (current.size>FFTNode.MIN_SIZE){
					numOfChildren=current.createChildNodes(childNodes);
					for (int i=0; i<numOfChildren; i++){
						stack[top++]=childNodes[i];
					}
				}
			}
		}
		/*
		for (int i=0; i<stack[0].cwm.length; i++){
			double error=Math.abs(real[i]-stack[0].cwm[i].re)+Math.abs(image[i]-stack[0].cwm[i].im);
			if (error>1.0e-4){
				throw new RuntimeException("error:"+error); 
			}
			real[i]=(float)stack[0].cwm[i].re;
			image[i]=(float)stack[0].cwm[i].im; 
		} 
		*/
		return stack[0]; 
		//System.out.println(f.isAllFree()); 
		//System.out.println("===>"+count); 
	}
	public static void testCorrectness(){
		/*
		for (int i=4; i<10; i++){
			check(i); 
		}
		*/
		//
		check(20); 
		 
	}
	public static void check(int order){ 
		int n=1<<order; 		
		System.out.println("size:"+n); 
		DepthFirstFFT fft=new DepthFirstFFT(n,order,16);
		boolean[] x=new boolean[n];
		float[] real=new float[n]; 
		float[] image=new float[n];
    	Random r=new Random(System.currentTimeMillis()); 
    	for (int i=0; i<n; i++){
    		x[i]=r.nextBoolean();
    		//System.out.println((x[i]?1:0)); 
    	}
    	long t0, t1; 
    	t0=System.currentTimeMillis(); 
		fft.fft(x,real,image);
		t1=System.currentTimeMillis(); 
		System.out.println("fft:"+(t1-t0)); 
		
        double[] fx=new double[n*2];
        for (int i=0; i<n; i++){
        	fx[i]=x[i]?1:0; 
        }
        
        if (false){
        	Complex[] c = new Complex[n];

        	// original data
        	for (int i = 0; i < n; i++) {
        		c[i] = new Complex(x[i]?1:0, 0);
        	}
        
        	c=FFT.fft(c);
        	for (int i=0; i<Math.min(n,16); i++){
        		double err=Math.abs(c[i].re-real[i])+Math.abs(c[i].im-image[i]);
        		if (err>1.0e-3){
        			for (int k=0; k<n; k++){
        				System.out.println(real[k]+"==>"+image[k]); 
        			}
        			FFT.show(c, "c"); 
        			throw new RuntimeException("err:"+err+" when n="+n +" at line "+i); 
        		}
        	}
        }
        System.out.println("OK"); 
        //FFT.show(c, "c"); 
        
	}
	public static void main(String[] args){		
		//testperformance(); 
		testCorrectness(); 
	}
	public static void testperformance(){ 
		int order=10; 
		int n=1<<order; 
		DepthFirstFFT fft=new DepthFirstFFT(n,order,16);
		boolean[] x=new boolean[n];
		long t; 
		long t0;

		int it=10000; 
		float[] real=new float[n]; 
		float[] image=new float[n];
		
		/*
		for (int i=0; i<it; i++){
			fft.fft2(x,real,image); 
		}
		
        long te=System.currentTimeMillis();
        System.out.println("fft2:"+(te-t));
		*/
        float[] fx=new float[n*2];
        for (int i=0; i<n; i++){
        	fx[i]=x[i]?1:0; 
        }
		
        FloatFFT_1D fft1D=new FloatFFT_1D(n);
        t=System.currentTimeMillis();
        for (int i=0; i<it; i++){
        	fft1D.realForwardFull(fx); 
        }
        
        long t3=System.currentTimeMillis();
        System.out.println("benchmark:"+(t3-t));

		t0=System.currentTimeMillis();
		for (int i=0; i<it; i++){
			fft.fft(x,real,image); 
		}
		t=System.currentTimeMillis();
		System.out.println("fft :"+(t-t0));

		/*
		long t; 
		long t0=System.currentTimeMillis();
		for (int i=0; i<10000; i++){
			FFTNode[] n=new FFTNode[1024*4]; 
			//FFTNode root=new FFTNode(0,1,8);
			for (int k=0; k<n.length; k++){
				n[k]=new FFTNode(0,1,8);
			}
		}
		t=System.currentTimeMillis();
		System.out.println("t :"+(t-t0));
		*/
	}
}
